(NASA-CR- 120889) STABILITY OF COMBUSTORS N72-20768 

WITH PARTIAL LB NOTH ACOUSTIC LINERS C.E. 

Mitchell, et al (Colorado State CJniv.) 

Mar. 1972 73 p CSCL 21H Unclas 

G3/28 23189 


\ 


\ 



NASA CR-120889 


ABILITY of COMBUSTORS tn 

length acousSc lS 


PARTIAL 


c - E. Mitchell, 


by 

W ' R ‘ ^spender and M. R. 


Baer 


COLORADO S1A TE UNIVERSITY 


Prepared for 

national aeronautics and space AD„t 

pace administration 


NASA Lewis p QO 

18 Resea «h Center 

Grant 06-002-095 
R1Chard J ' ^eet Ma„ ager 


NOTICE 


This report was prepared as an account of Government- 
sponsored work. Neither the United States, nor the 
National Aeronautics and Space Administration (NASA), 
nor any person acting oiji behalf of NASA: 

A. ) Makes any warranty pr representation, 

expressed or implied, with respect to the 
accuracy, completeness, or usefulness of the 
information contained in this report, or that 
the use of any information, apparatus, method, 
or process disclosed in this report may not 
infringe privately-owned right; or 

B. ) Assumes any liabilities with respect to the 

use of, or for damages resulting from the 
use of, any information, apparatus, method 
or process disclosed in this report. 

As used above, "person acting on behalf of NASA" 
includes any employee or contractor of NASA, or 
employee of such contractor, to the extent that such 
employee or contractor of NASA or employee of such 
contractor prepares, disseminates, or provides access 
to any information pursuant to his employment or 
contract with NASA, or his employment with such 
contractor. 


Requests for copies of this report should be referred to 

National Aeronautics and Space Administration 
Scientific and Technical Information Facility 
P. 0. Box 33 

College Park, Md. 20740 



1. Repo 1 1 No. 2 Government Accession f Jo. 

NASA CR 120889 


A, Title and Subtitle 

Stability of Combustors with Partial Length 
Acoustic Liners 


3. Recipient's Catalog No. 


G. Reix>rt Date 

March 1972 


6. Performing Organization Code 


7. Authors) 

C. E. Mitchell, W. R. Espander and M. R. Baer 


8. Performing Organization Report No. 


>0. Work Unit No. 

9. Performing Organization Name and Address 

Colorado State University — — r — rrr- 

J 11. Contract or Grant No. 

Fort Collins, Colorado 80521 N( ^ R Qg_oo2-095 

13. Type of Report and Period Covered 

12. Sponsoring Agency Name and Address Contractor Report 

National Aeronautics and Space Administration 14 . Sponsoring Agency Code 
Washington, D. C. 20546 

15. Supplementary Notes 

Project Manager, Richard J. Priem, Chemical Propulsion Division, 

NASA Lewis Research Center, Cleveland, Ohio 


16. Abstract 

An analytical technique for the evaluation of combustion stability 
in rocket motors with partial length acoustic absorbers is presented. 
The combustors considered in this work have concentrated combustion 
zones at the injector, finite mean flows, cylindrical cross sections, 
and acoustic liners of arbitrary length and impedance. Linear three 
dimensional oscillations in such combustion chambers are analyzed 
using an integral equation-iteration technique. Stability limits in 
terms of a combustion response factor are calculated for several 
values of Mach Number, length to radius ratio, liner impedance, 
liner length, liner location, and nozzle admittance. Results indi- 
cate that increasing liner length increases combustor stability sub- 
stantially at low Mach Numbers but has a substantially smaller effect 
for larger Mach Numbers. Increasing Mach Numbers or length to radius 
ratio have destabilizing effects while liner location has only a 
minor effect on stability. 


17. Key Words (Suggested by Author (s)) 

Combustion instability 
Liquid rockets 
Acoustic absorbers 


18. Distribution Statement 


Unclassified - unlimited 


19. Security Classif. (of this report) 

Unclassified 


20. Security Oassd. (of this page) 

Unclassified 


21. No. of Pages 

71 


For sale by the National Technical Information Service. Spnngfield. Virginia 22151 

If 











preceding page blank 


iii 


Foreword 

The research described herein, which was conducted at 
Colorado State University during the period December 1, 1969 
to December 31, 1971 was supported by NASA Grant NGR 06-002-095. 
The work was done under the management of the NASA Project 
Manager, Dr. Richard J. Priem, Chemical Rockets Division, NASA 


Lewis Research Center. 



ABSTRACT 


An analytical technique for the evaluation of combustion 
stability in rocket motors with partial length acoustic absorbers 
is presented. The combustors considered in this work have con- 
centrated combustion zones at the injector, finite mean flows, 
cylindrical cross sections, and acoustic liners of arbitrary 
length and impedance. Linear three dimensional oscillations in 
such combustion chambers are analyzed using an integral equation- 
iteration technique. Stability limits in terms of a combustion 
response factor are calculated for several values of Mach Number, 
length to radius ratio, liner impedance, liner length, liner 
location, and nozzle admittance. Results indicate that increas- 
ing liner length increases combustor stability substantially at 
low Mach Numbers but has a substantially smaller effect for larger 
Mach Numbers. Increasing Mach Number or length to radius ratio 
have destabilizing effects while liner location has only a minor 
effect on stability. 
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Nomenclature 

speed of sound 

- quantities defined on Page A-2 

- quantity defined on Page A-2 

- functions defined on Page A-4 

- Green’s function 

- modified Green’s function defined on Page A-2 

- unit imaginary = /-T 
imaginary part of ( ) 

- positive integer 

- Bessel Function of the first kind of order m 

- absorber impedance 

- positive integer 
chamber length 

- Mach Number 
positive integer 

- positive integer or interaction index 

- outward unit normal vector 

- pressure 

- velocity vector 
chamber radius 

- radial length 
position vector 
source position vector 

- real part of ( ) 
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- acoustic admittance 

- ratio of specific heats 
acoustic eigenvalue for mode N 
circumferential coordinate 

- imaginary part of complex frequency 

root of J (X 0 ) = 0 
m £m 

- normalizing constant defined on Page A-l 

matrix discussed on Page A-3 
density 

sensitive time lag 

- perturbation velocity potential 

- normalizing constant defined on Page A-2 
acoustic eigenfunction for mode N 
complex frequency of oscillation 
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Subscripts 

C - quantity evaluated on chamber cylindrical periphery 

i - injector quantity 

N - nozzle quantity 

o - mean value 

Superscripts 

- perturbation quantity or derivative with respect to 
argument 

* - dimensional quantity 

- - chamber quantity with no absorber present 

- designates particular mode integers chosen 

(j) - approximation of order j 
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INTRODUCTION AND SUMMARY 


Acoustic liners designed for liquid rocket combustors usually 
cover only a fraction of the available wall area. It is therefore de- 
sirable to have an analytical means of predicting the stability behavior 
of combustors with acoustic absorbers of this type in addition to the 
existing analyses which are restricted to fully lined or unlined chambers 
It is the purpose of this report to present an analytical technique of 
this type and to demonstrate its efficacy in predicting the stability 
characteristics of combustors with partial length liners. Priem (1) has 
shown that the combustion zone response, the nozzle admittance, and the 
mean flow in the chamber all can have strong effects on the frequency 
of chamber oscillations and on the linear stability of the oscillations 
both for unlined and fully lined chambers. In addition he showed that 
the axial dependence of the wave amplitude and phase was quite sensitive 
to the parameters mentioned above. Consequently, it is clear that in 
order to get a complete picture of the stability behavior of a combustor 
with a partial length liner, it is necessary to include mean flow, nozzle 
and combustion zone response effects in any partial length liner analysis 
Because of the discontinuous boundary condition introduced into the 
problem by the partial length liner, standard separation of variables 
techniques will not work in the solution of the partial length liner 
problem. Instead it is necessary to transform the partial differential 
equations governing the chamber flow field into integral equations using 
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Green’s theorem and a Green’s function for the chamber. Similar trans- 
formations have been used previously by Culick (2) and Oberg (3) . 

These integral equations are then solved using an iterative technique 
and yield results in terms of stability limits and pressure and veloc- 
ity waveforms. The effects of mean chamber Mach Number, nozzle imped- 
ance, liner impedance, length and location, chamber length to radius 
ratio and combustion zone response are then investigated using this 
iterative solution method. 
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THEORY 

Combustion Chamber Model 

The type of combustor considered in this work is assumed to have a 
concentrated combustion zone at the injector, to be of cylindrical cross 
section, to be terminated by a supercritical nozzle of known impedance 
and to have a finite mean flow (mean chamber Mach Number) . An arbitrary 
fraction of the cylindrical periphery is assumed to consist of an acoustic 
liner of known impedance. The rest of the cylindrical chamber surface is 
assumed to have an infinite impedance. Figure (1) shows pictorially the 
type of combustor considered. 

The response of the combustion zone mass generation rate is assumed 
to be pressure dependent only. Therefore, it can be represented using 
either a combustion response factor of the type used by Priem (1) or by 
the sensitive time lag (n, x) model of Crocco (4). Both models will be 
used in interpreting the results of calculations performed in this work. 
Downstream of the combustion zone the gasdynamic field is assumed to be 
composed of a single component, calorically perfect gas. The flow is 
taken to be homentropic and irrotational . 

Because of the irrotationality of the flow it is possible to com- 
bine the conservation equations into a single partial differential equa- 
tion in terms of the velocity potential, $>, following standard techniques. 
(See for example Ref. 1). If the oscillations are considered to be of 
small amplitide and periodic in time then $ = Mz + (jje 10 ^, where M is the 
mean Mach Number for the flow, z is the axial coordinate nondimension- 
alized by the chamber radius, cf) is the perturbation velocity potential, 
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i is the square root of -1, 03 = to + iX is the complex frequency of 
the oscillations and t is time nondimens ionalized by the ratio of the 
chamber radius to the mean sound speed. With these definitions the 
linearized partial differential equation for <J> can be written 

v 2 cj> + 0 ) 2 (j> = M 2 |— |- + 2i03 yk (1) 

dZ z a z 

Perturbation pressure and velocity are related to <f> by the following 
expressions 

p = - yio3<j) - yM > q = V<j) 
dz 

The boundary condition on Equation (1) is V<J> • n = Sp on the surface of 

the chamber (including injector plane and nozzle exit plane). 3 is the 

specific acoustic admittance of the bounding surface in question. On 

the portion of the cylindrical walls formed by the acoustic liner 

3=3 = 1/yK, where K is the specific acoustic impedance of ‘the liner, 

c 

At the combustion zone (injector plane) 3 = 3^ = M(— - - N) , where N is 

m 1 

the combustion response factor defined as N = — and y is the specific 
heat ratio, m is the normalized perturbation in mass flux. If the 
time lag model is used 3^ = M[n(l-e^ WT ) - y] where n is the interaction 
index and T is the time lag. At the nozzle entrance plane 3 = 3^, the 
nozzle admittance. In this work the nozzle admittance used was either 
that for a short nozzle or the value obtained using the tabulated results 
of Crocco and Sirignano (5) for conical nozzles. 
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Basic Integral Equations 

The single partial differential equation for the chamber flow 
(Equation (1)) and the boundary conditions discussed above can be trans- 
formed into two integral equations through the use of Green’s theorem 
and the introduction of a Green’s function for the chamber. The result- 
ing equations are 


<f> = + f + YM If] ds o 

■'S*' 


'///• 


+ 2Mia) ^ dV 


dz O 


2 2 
or - n 


■fffi 




M » f t~t + 2MifcO dV 

N / / / N dz Z dz O 


//* 


+ / /Bfyyiwct) + YM f|) ds o 


( 2 ) 


(3) 


In Equations (2) and (3) the quantities Tl^ and G^ have the following 

meanings, is one of the acoustic eigenfunctions for the chamber with 

no mean flow, combustion zone, liner, or nozzle. T)^ is the eigenvalue 

for the frequency of oscillation associated with the eigenfunction Q . 

N 

G^ is a Green’s function for the chamber represented as an eigenfunction 
expansion in the normal acoustic eigenfunctions of the chamber, with one 
specific eigenfunction deleted from the triply infinite sum. The 

exact forms of the eigenfunctions, eigenvalues and Green’s function are 


given in Appendix A. 
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Solution of the Equations 

It was desired to obtain results in the form of stability limits 


in terms of either the real and imaginary parts of the combustion re- 
sponse factor of in terms of n and T. In order to obtain results of 
this type Equation (3) was rewritten in the following form where g^, 
the complex combustion admittance appears explicitly on the left hand 


side of the equation. 


2 2 
or - n 

'n 


-f/M> + 2Miw H 5 dv 0 ■ + ym H ) ds o 

V s-s-j 

ffa* + YM ds o 


In this equation s^ is the injector plane surface area and s-s^ is the 
rest of the chamber surface area. 


Equations (2) and (4) above are sufficient to determine <j> and 
for given mean chamber operating conditions and geometry, and liner length, 
impedance, and location for a particular value of the complex frequency, 
to. In order to obtain stability limits in terms of g_^, the equations 
must be solved for a range of frequency values. Stability limits in 
terms of the combustion response factor or n and T are easily obtained 
using the expressions relating these quantities to g given earlier. 

For a given frequency the following iterative solution method was 
followed. First, an initial guess for the form of <J> and the value of g^ 
which will satisfy Equations (2) and (4) was made. In this work the <j) 
and g^ resulting from solution of Equation (1) for the case of an un- 
lined combustor were used as the initial approximations. This solution 
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is easily obtained using the technique of separation of variables since 
no discontinuous boundary condition is present in this case. These 
trial values for cf> and 3^ were then substituted into the integrals on 
the right hand sides of Equations (2) and (4), resulting in first approxi- 
mations for <f> and 3^ for the lined chamber configuration. These first 
approximations were then themselves substituted into the right hand side 
integrals to produce improved approximations and so on. The iteration 
process continued until successive approximations were suitably invarient. 
A one per cent difference in successive values of 3^ was then taken to 
be the cut off point in most of the calculations performed, though this 
degree of accuracy is, of course, arbitrary. 
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RESULTS AND DISCUSSION 

Calculations were performed for several values of Mach Number, 
chamber length to radius ratio, liner length, liner impedance, liner 
location and nozzle admittance. Results of those calculations in terras 
of both neutral stability limits and pressure and velocity profiles 
will be presented. Results for the first transverse mode of oscilla- 
tion only will be discussed although some second and third transverse 
calculations have been done with similar results. Before giving these 
results a brief discussion of the convergence and accuracy of the 
iterative technique used will be presented. 

It was found that the iteration technique converged in less than 

/ 

ten iterations to within the specified accuracy for almost all combina- 
tions of chamber design parameters investigated with two exceptions. 

The first of these was related to the magnitude of the damping effect 
of the liner. If this damping effect was too large (typically of the 
order of 100% decay per cycle or 1500 sec 1 damping coefficient) con- 
vergence was poor or did not occur at all. A rule of thumb seemed to 
be that if the ratio of liner length to liner impedance was less than 
about 0.50 convergence was good. Above this value difficulties began 
to occur. For a one third length liner, for example, the calculations 
were restricted to liners with impedances in excess of 0.60. This is 
clearly a very small impedance for practical liners so the restriction 
does not seem important for realistic calculations. 

The other situation in which convergence was a problem was for the 
part of the neutral stability curves occurring in the transition between 
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various pure and combined modes of oscillation. That is, for example, 
for pure first transverse mode oscillations or for combined first trans- 
verse, first longitudinal oscillations convergence was good. For fre- 
quencies corresponding to oscillations not clearly first transverse or 
first transverse, first longitudinal, difficulties in convergence were 
experienced. These transition regions appear as loops on the neutral 
stability curves. Points along the loops can be found by improving 
the initial guess for cf> and 3^ with the inclusion of both modes in 
question, though convergence is still slower than usual. 

Because the Green’s function, G^, is represented as an eigenfunction 
expansion, <|)^\ the jth approximation for the velocity potential is also 
represented as a doubly infinite series in the r and z directions. (See 
Appendix A) . Explicitly 


,«> - i + 


OO 00 

<t> + cos m9 y, y] 

n=l 


(i) T \ nTtz 

u ln J " (X Jm r) cos T 


where (p is the analytical solution for cp when no liner is present, m 
signifies the particular transverse primary mode being considered (e.g., 

/v s (i) 

m = 1 denotes the first transverse mode) and y„ is the coefficient 

Jen 

matrix for the doubly infinite series. In order to carry out numerical 
calculations must be finite in dimension. The largest matrix used 

in this work had 10 terms in the £ series (determining behavior in the 
radial direction) and 100 terms in the n series (determining behavior 
in the axial direction. Increasing matrix size beyond 10 x 100 can be 
shown to affect the stability limit results by less than 3%. Table 1 


shows the errors in combustion response factor incurred when smaller 
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matrices are used as well as the number of iterations needed for con- 
vergence and the computer time required (CDC 6400). For the calcula- 
tions performed in this study a 3 x 30 matrix was used in most cases. 
The effect of matrix size on pressure and velocity profiles will be 
discussed later. 

Stability Limit Results 

Convergence of the technique in terms of shape and location of the 
linear stability curve is shown in Figure (2). Re(N) and Im(N) are the 
real and imaginary parts of the combustion response factor, M is the 
mean Mach Number, L/R is the length to radius ratio for the chamber, K 
is the acoustic impedance of the liner, and G is the acoustic response 
of the nozzle. A value for G of 0.9166 corresponds to a short nozzle 
response. The normalized liner length is given by X. A value of X of 
2.7 indicates a fully lined chamber if L/R = 2.7. The nondimens ional 
frequency co /to q is the parameter used along the curves. u) is the actual 
chamber frequency; is the frequency of the acoustic first transverse 
mode for no combustion, mean flow, or liner. Progressing along any of 
the curves in the direction of increasing frequency, regions to the 
right of the curve are linearly unstable while regions to the left are 
linearly stable. For the example chosen it is seen that there is no 
change in the stability limit location after five iterations. It is 
to be noted that the first and second approximations are noticeably in 
error even though M = 0.10 and the usual linear or second order expan- 
sions in Mach Number might be expected to be accurate. 
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Figure (3) shows the effect of varying liner length on the stability 
of a low mean Mach Number chamber. It is readily seen that increasing 
the liner length increases the region of linear stability. (The neutral 
stability curves shift to the right with increasing X). Moreover, this 
increase in stability with liner length is seen to be almost linear for 
the Mach Number being considered. In Figure (4) which is the same type of 
plot, the effect of increasing the mean Mach Number is seen. For this 
figure all chamber parameters are the same with the exception of the mean 
Mach Number which has been increased from 0.1 to 0.33. It can be seen that 
"loops" are introduced into the neutral stability curves. These coincide 
with the appearance of combined modes of oscillation and have been observed 
previously by Priem (1) and others. The stabilizing effect of adding a 
liner is diminished with the increase in Mach Number. As can be seen, the 
shift to the right caused by the addition of a liner is considerably less 
than it was for the case of the low Mach Number chamber in Figure (3) . 
Finally it is seen from Figure (4) that adding more liner is also less 
effective than was the case for M = 0.1. Indeed, a 1/ 3 length liner pro- 
vides half the stability improvement of a full liner and a 2/3 length 
line provides essentially the same stabilizing effect as a full liner. 
Figure (5) shows the effect of decreasing the impedance (increasing the 
damping effect) of the liner while keeping all other parameters at the 
values they had in Figure (4). Two facts are apparent from the figure. 
First the linear stability is improved considerably as all three liner 
neutral stability curves are shifted to the right from their position 
in Figure (4). Second, the .effect of liner length has become even less 
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important. Indeed, there seems to be little difference between a 1/3, 
2/3 and full liner as far as overall stability improvement. These con- 
clusions are verified and perhaps clarified in Figure (6) . This figure 
shows, for the same chamber parameters and the same liners, neutral 
stability plots on the n,T plane, n is the interaction index and T 
the sensitive time lag of Crocco's Sensitive Time Lag Model (4). Fre- 
quency increases from the right to the left along the curves. The 
region above any neutral stability curve is a region of linear stabil- 
ity, the region below a region of linear stability. Thus, shifting a 
neutral stability curve upward increases stability. Figure (6) shows 
that for lower frequencies the longer liners are slightly more stabiliz- 
ing, but that for higher frequencies (e.g., combined 1st transverse 1st 
longitudinal) the shorter liner is actually slightly more effective. 
Figures (7) and (8) show, respectively, the effect of increasing Mach 
Number for a fixed liner impedance and length and the effect of decreas- 
ing the liner impedance for fixed Mach Number and liner length. It is 
clear from the figures that increasing the Mach Number is destabilizing 
while decreasing the liner impedance provides a stabilizing effect. 

Figure (9) shows the effect of increasing the chamber length to 
radius ratio on the chamber’s stability. It can be seen that increasing 
L/R provides a substantial destabilizing effect. Examination of Figures 
(7) and (8) indicate that increasing either L/R or M increases the num- 
ber of combined modes possible. Figure (10) demonstrates the influence 
of liner location. A 1/4 length liner is located at various positions 
in the chamber and the neutral stability curve determined for each 
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location. If Xi = 0 the liner is at the injector end, if Xi = 0.75 the 
liner is at the nozzle end, and if Xi = .40 the liner is in the center 
of the chamber. It can- be seen that there is little difference in over- 
all effectiveness between Xi = 0.0 and Xi ® 0.40 but that Xi = 0.75 is 
somewhat less effective. The influence of nozzle response on the neutral 
stability limits is shown in Figure (11) . Two types of nozzles are com- 
pared both with and without 1/3 length liners. One of the nozzles is a 
"short” nozzle, the other is a conical nozzle of the type analysed by 
Crocco and Sirignano. For the short nozzle the nozzle admittance is not 
a function of frequency and provides a small damping effect. For the 

conical nozzle is a function of a), and, according to the Crocco- 
N / 

Sirignano calculations, provides a small driving effect in the frequency 
range of interest. Comparison of the short nozzle and conical nozzle 
stability limits for either lined or unlined chambers indicates that, 
as is to be expected, the chamber with the short nozzle is always more 
stable than the chamber with the conical nozzle. It is important to 
observe that the stability improvement caused by the addition of a liner 
is approximately the same for both types of nozzles. This means that 
stability improvement predictions for lined chambers with short nozzles 
are likely to also apply for other types of nozzles. 

Pressure and Velocity Profiles 

The perturbation pressure and velocity were calculated as a function 
of spatial coordinates using the relationships given in Appendix A. Re- 
sults are presented for a frequency value (oo/oo^) of 0.90. Figure (12) 
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is a plot of relative pressure amplitude (ratio of pressure amplitude 
at a particular radial and axial location to the pressure amplitude 
at the injector for the same radial location) as a function of axial 
length for the standard chamber with no liner , with a 1/3 length liner 
(X = 0o9) and with a full liner (X = 2.7). 

For the unlined chamber the amplitude continuously decreases from 
injector to nozzle. For the fully lined chamber the amplitude decreases 
very rapidly near the injector and then maintains a more or less constant 
value for the rest of the chamber. The 1/3 length liner causes a more 
rapid decrease in amplitude near the injector than in the case of the 
unlined chamber, but not as pronounced as in the case of the fully lined 
chamber. Approximately at the end of the lined portion the 1/3 length 
liner curve starts to follow the unlined chamber curve and converges to 
it as the nozzle end is approached. 

Figure (13) shows the phase angle difference between the injector 
and a given location for a given wave. The same chamber configuration 
and liners are considered as in Figure (12). It can be seen that for 
the unlined chamber the wave downstream of the injector always leads 
the wave at the injector with the lead in general increasing with axial 
length. For the chamber with a 1/3 length liner the phase angle curve 
is qualitatively similar except near the injector where a region of 
lagging occurs. The fully lined chamber presents a quite different 
wave shape with the wave in the center of the chamber lagging the wave 
at both the injector and nozzle. The effect of matrix size on the 

relative pressure amplitude curve is shown in Figure (14). The large 
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matrix is 10 x 30 and the small matrix is 3 x 30. The curves are 
qualitatively very similar though the small matrix predicts a somewhat 
more rapid decrease in pressure amplitude than does the large matrix. 

The radial velocity as a function of axial length at a radial 
location 0.90 of the chamber radius is shown in Figure (15) for the 
same chamber with the same three liner configurations. Velocity values 
were calculated assuming that the arbitrary amplitude factor for the 
velocity potential* (j), was unity in the case where no liner was present. 
It can be seen that, as might be expected, (since v^ at the wall 4 0) 
the addition of a full length liner increases the radial velocity above 
the value with no liner present at every axial location. The 1/3 length 
liner causes an even greater increase in v^ over the lined portion but 
v^ gradually tends toward v^ for the unlined chamber as the nozzle is 
approached. 

The dependence of the radial velocity on radius for an axial loca- 
tion 0.45 of the chamber radius is shown in Figure (16). As the wall 
is approached (r = 1) the radial velocity increases substantially above 
the value for the unlined chamber and then decreases rapidly to zero 
both for the 1/3 length liner and the full liner. The 1/3 length liner 
produces slightly larger radial velocities than does the full length 
liner. Because of the use of the Green’s Function Technique, the wall 
boundary condition v^ 4 0 cannot be satisfied exactly at the wall ex- 
cept in the limit as an infinite number of terms are retained in the 
radial direction in the matrix for y^^ . 

A comparison of small matrix (3 x 30) to large matrix (10 x 100) 
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results near the wall is shown in Figure (17). As the number of terms 
kept is increased a better satisfaction of the wall boundary condition 
is experienced. It is to be recalled that the increase in accuracy in 
the determination of combustion response factors improved by less than 
5% in going from the 3 x 30 to 10 x 100 matrix. Therefore, it seems 
that the detailed satisfaction of the wall boundary condition is rela- 
tively unimportant in obtaining neutral stability results. Figure (18) 
is similar to Figure (16) except that the axial location has been 
changed to 1.8 of the radius. Here the wall is unlined for the 1/3 
length liner chamber. It can be seen that for this case the fully lined 
chamber radial velocity is larger than the other two near the wall and 
that the curve for the 1/3 lined chamber is considerably closer to the 
unlined chamber curve than it was in Figure (16) . 
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CONCLUSIONS 

An analytical method for predicting the linear stability behavior 
of rocket combustors with partial length liners has been presented. 

It has been shown that the method is applicable over a wide range of 
chamber and absorber design parameters. Calculations performed indi- 
cated that the following general characteristics of the stability of 
combustors with partial length liners are to be expected. 

1. An increase in the mean flow Mach Number provides a de- 
stabilizing effect for a given chamber and reduces the 
stability improvement provided by an acoustic liner. 

2. Increasing the chamber length to radius ratio results in 
the same type of destabilizing effects as increasing the 
Mach Number. 

3. Increasing liner length as a means of increasing stability 
is most effective at high liner impedances and low Mach 
Numbers. For large Mach Numbers and small impedances 
partial length liners can be as effective as full length 
liners . 

4. Decreasing liner impedance for a given liner length always 
provides a marked increase in stability. 

5. Location of a partial length liner is not very important, 
though liners located at the nozzle tend to be slightly 
less effective than if located elsewhere. 

6. A wide range of frequencies of oscillation (including 
combined modes) is possible with a partial length liner 
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in place. The liner is more effective as a damping 
device (even if its impedance is fixed) at some fre- 
quencies than it is at others. 
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Figure 7 Effect of Mach Number on neutral stability - 1/3 length liner 
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Figure 12 Pressure profiles for various length acoustic liners. 
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Figure 14 Comparison of small matrix and large matrix pressure profiles. 
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Figure 15 Radial velocity profiles for various length acoustic liners. 
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Figure 16 Radial velocity profiles for various length acoustic liners. 
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Figure 18 Radial velocity profile for various length acoustic liners, Z = 1.80 
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APPENDIX A 
Basic Equations 


The partial differential equation describing the chamber oscilla- 
tions is 

V 2 (p + o) 2 c p = M 2 + 2Mia) — 
r r dz 2 dz 

with the boundary condition V<() • n = - ySCioje)) + M ^-) on all surfaces. 

d z 

In order to convert this equation into an integral equation a Green’s 
Function which satisfies the following equation is introduced 


V 2 G + o) 2 G = 6 (r-r Q ) 

where VG • n = 0 on all boundaries and 6 is the usual delta function. 
The Green’s Function is represented as an expansion in the orthonormal 
eigenfunctions for the chamber with no combustion or mean flow and solid 
walls as 


G = 


fi £mn (r o ) 


o co 2 - rio 

x, m n Jcmn 


where is one of the orthonormal eigenfunctions given by 


and 


fL 

£mn 


cos m0 J (X n r) 
m £ro 


cos 


nTT z 
L ■ 


A 


_ 1 _ 

Zmn 


\mn 


n 2 TT 2 

L 2 


+ X 


2 

Z m 


The quantity X 2 is a root of the equation J (X n r) » 0 
Zm m Zm 

L Jr=l 


and A 


£mn 


is determined by the condition 


///« 


ft* dV = 1 
Zmn 


The modified Green’s Function, G^, which appears in the integral 
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Equations (3)aand (4) is just the Green's Function with one term in the 
series removed. That is, the term with SL = £ , m = m, and n = n (where 

Thus 


A A 


iy m and n are arbitrary) is deleted from the series for G^. 


(r) (r ) 

„ _ , V ^&mn ^£mn (1 - 6 (l ,£)6 (m,m)6 (n,n)) 

G „ = — -_" _ n2 


N 


r v 

H m n 


CO 


£mn 


where 6(i,j) = 1 if i = j and fi(i»j) * 0 if i ^ j, 

In Equation (2) the quantity can be interpreted as and in 

Equation (3) rj.. is equivalent to In the solution of Equations 

N x,mn 

(2) and (3) the velocity potential and combustion response for the cham- 
ber with no liner but with mean flow are used as first approximations 
for cj) and g^. These quantities are called respectively $ and g^. The 
form of cj) can be found using the technique of separation of variables 
and is given by 


where' 


$ 


Bi 


B 2 


cos 


fM 

£mn 


me T ,, v f izBj 

— J s ( ha r)(e 


+ c 


e lzBj ) 


<P 


cjM + [o) 2 M 2 + (M 2 -1)(X|/v - w 2 )]^ 
(1-M 2 ) 

toM - [u) 2 M 2 + (M 2 -1)(X|~ - w 2 )]^ 
_____ 


giLBi [Bj + 6 n ( Y o) + yMBi )] 
e iLB 2 [B 2 + $ N (yw + yMB 2 ] 


and ij; is determined by the condition 
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then determined from the boundary condition at the injector, 
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or 


*L 

3z 

(Yi<W> + YM ||) 


Substitution of $ and 8^ into Equations (2) and (3) as the initial step 
in the iteration procedure described earlier eventually leads to the 
following two relationships for the calculation of the j*-* 1 approxima- 
tions for 4> and 8^ 


,( 3 > .2 


= * + cos S J m (A £m r) C0S 


SL n 




8 


(j+l) ... 


N 


/* b 
y n 


f 2 (<f> (j) )ds 


S . 
1 


where is a matrix dependent only upon <J>^ and 8^"^ (y^^ depends 

on $ and 8^» 8^ depends on $ and 8.), f 2 (y) = Yiwy + M 

IX 1 oz 

3 v 3 ^ v 

f i (y) = 2iw + M - 5 — 5 -, s. is the surface of the injector, s is the 

0 z 0 Z 1 IN 

surface of the nozzle exit plane, is the surface of the liner and V 

is the chamber volume. For any solution, values of m, and n must be 

selected. These integers will determine the general acoustic waveform 

of <J) expected and the exact waveform for <p . For example m = 1, * = 1, 

n = 0 gives the solution for the pure first transverse mode. If the 

AAA 

wrong values of m, Jo, n are selected for a given frequency value the 
iteration technique will fail. That is, if the solution is a first 
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transverse mode type oscillation and a combined first transverse - second 
longitudinal mode is assumed (by taking m = 1, SL = 1, n = 2) the itera- 

A A 

tion scheme will not work until the value of n is changed to n = 0. 
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APPENDIX B 
Computer Program 

This section describes the use of the computer program, COMBADM, 
for the determination of the neutral stability curves for a given set 
of input parameters which characterize the combustion chamber. The 
program was written in FORTRAN IV to run on a CDC 6400. Very little 
effort was made to optimize the program with respect to running time. 

The program, COMBADM, calculates two neutral stability curves 
for a given combustor as a function of frequency. The first combustion 
admittance calculated for each frequency is the neutral stability point 
for the combustor with no liner. This point corresponds to BETAWI in 
the output. The last value of the combustion admittance listed for 
each frequency is for the combustor with the acoustic liner. This 
point is determined by the method of successive approximations. The 
first approximation to the combustion admittance for the lined chamber 
is that of the unlined combustor. The combustion admittance is then 
used to calculate a velocity potential for the lined combustor and 
then a new admittance is calculated for the combustor. This technique 
proceeds until the error between successive iterations is within pre- 
determined limits. The accuracy of the combustion admittance calculated 
is limited only by the size of the coefficient matrix used in the cal- 
culation of the velocity potential. 

The input parameters are: 

1. LENGNO - the problem identifier 

2. F - the initial frequency ratio of the combustor 

3. ALENGTH - the length to radius ratio for the combustor 
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4. AMACH 

5. GAMMA 

6. G 

7. AK 

8. XI 

9. X2 

10 . ERROR1 

11 . ERROR 2 

12 . LHAT 

13. MHAT 

14 . NHAT 

15 . LDEX 

16 . NDEX 

17 . MATRIX 

18. FINC 


- mean flow Mach Number 

- ratio of specific heats 

- complex number specifying the acoustic impedance 
of the nozzle 

- complex number specifying the acoustic impedance 
of the liner 

- the beginning of the liner opening with respect 
to the injector 

- the end of the liner opening with respect to 
the injector 

- the maximum acceptable error in the real part of 
the combustion admittance 

- the maximum acceptable error in the imaginary part 
of the combustion admittance 

- radial acoustic mode number 

- transverse acoustic mode number 

- longitudinal acoustic mode number 

- specifies the number of radial terms in the 
coefficient matrix 

- specifies the number of longitudinal terms in 
the coefficient matrix 

- identifier used to determine whether or not the 
coefficient matrix is printed out 

- the increment size of the frequency ratio 
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19. FMAX - the maximum frequency for which a combustion 
admittance is calculated 

The output from the program lists the input variables, the com- 
bustion admittance and the combustion impedance. 

Sample Calculation 

An example of the output from COMBADM is shown on the page imme- 
diately following the listing of the program. The example shown is 
typical of many cases that can be run. In the output the neutral sta- 
bility point for the unlined combustor is denoted as BETAWI and the 
neutral stability point for the lined combustor to the specified accuracy 
is the last two numerical values listed, where the first of these values, 
BETAI(J), is the combustion admittance and the latter value, N, is the 
combustion impedance. The only problem likely to be encountered in the 
program is the selection of the appropriate longitudinal acoustic mode 
(NHAT) , i.e., as the frequency increases higher order longitudinal modes 
are more dominant in the solution and NHAT often must be increased from 
zero to one or even two in order for the iteration process to converge. 
(See Appendix A) 

The form of the input for COMBADM is shown immediately following 
the list of the program. The input is on three data cards as follows: 


Card 1 



Columns 

Variable 

Type 

1-10 

LENGNO 

Alphanumeric word 

11-20 

Re(F) 

Real number 


21-30 


Im(F) 


Real number 
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Columns 

Variable 

Type 

31-40 

ALENGTH 

Real number 

41-50 

AMACH 

Real number 

51-60 

GAMMA 

Real number 

61-70 

Re(G) 

Real number 

71-80 

Im(G) 

Real number 

Card 2 



Columns 

Variable 

Type 

1-10 

Re(AK) 

Real number 

11-20 

Im(AK) 

Real number 

21-30 

XI 

Real number 

31-40 

X2 

Real number 

41-50 

ERR0R1 

Real number 

51-60 

ERR0R2 

Real number 

61-63 

LHAT 

Integer 

64-66 

MHAT 

Integer 

67-69 

NHAT 

Integer 

70-72 

LDEX 

Integer 

73-75 

NDEX 

Integer 

76-80 

MATRIX* 

Alphanumeric word 


*If the coefficient matrix is to be printed out, then MATRIX is printed 
as YESBB, where B signifies a blank. 
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Card 3 


Columns 

Variable 

Type 

1-10 

Re(FINC) 

Real number 

11-20 

Im(FINC) 

Real number 

21-30 

Re(FMAX) 

Real number 

31-40 

Im(FMAX) 

Real number 

The last page is 

the output and includes 

the essential input 


variables and the two neutral stability points. 
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Nomenclature 



COMB ADM 


U) 

OMEGA 


M 

AMACH 


Y 

GAMMA 


n 

Function 

ETA 


Function 

BESPRIM 


Function 

BESSEL 

a L 

Function 

FNORM 

¥ 

Function 

PSI 

Bi 

BETA1 


b 2 

BETA2 


C 

C 


h 

BETAWI 



BETAI or 

TETAIS 


BETAC 


6 n 

BETAN 



AMU 
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PROGRAM COMBADMt OUTPUT* 101,TAPE6=OUTPUT, 1NPUT=101 ,TAPE5=I^UT) 
COMPLEX BETA 1 , BETA2, BETAN, BETAC, BETA I, BET AW I , OMEGA, C, AK1 , 

1 AK2, PS12, AMU, BETA I S , C12 

COMPLEX G, AN, AK, A, A4, A5, BETA, DUM1 , DUM2, PS!, DUMMY, Cl, 

1 G2,F,AKA,F]NC,FMAX 

COMMON / BLKA / BETA 1 , BETA2, AMACH, GAMMA, OMEGA, ALENGTH, C, AK! 
1, AK2, BETAN, RLMDA, BETA I, BETAC, BET AW 1 , LHAT, MHAT, NHAT, 

1 PS 1 2 , X 1 , X2/BLKB/AMU (10,30,29) .BETA IS (30) , A 
COMMON C 1 2 (30 , 30 ) , I SET , LDEX , NDEX , ERROR 1 , ERR0R2 

2 READ (5,200) LENGNO , F , ALENGTH , AMACH , GAMMA , G , AK , X 1 , X2 , ERROR 1 , ERR0R2 , 

1 LHAT , MHAT , NHAT , LDEX , NDEX , MATR I X , F I NC , FMAX 

IF ( EOF , 5 ) 50 , 3 

3 A = CMPLX ( 1.0 , 0.0 ) 

RLMDA = BESPRIM ( MHAT + 1 , LHAT ) 

7 CONTINUE 
NHAI = NHAT + 1 
OMEGA-F’ 1.84 11 8378 
BETAC = 1./GAMMA/AK 
BETAN = AMACH * ( G - 1. / GAMMA ) 

BETA=CSQRT (OMEGA* AMACH ‘OMEGA ‘AMACH+ (AMACH* AMACH- 1 . ) * (RLMDA* *2 - 
IOMEGA ’OMEGA) ) / ( 1 . -AMACH* AMACH) 

BETA1 = (OMEGA * AMACH) / ( 1 . -AMACH* AMACH) +BETA 

BETA2* (BET A 1 ) -2. ‘BETA 

A4 = CMPLX ( -A I MAG (BETA 1) , REAL (BETA 1 ) ) 

DUM1 = CMPLX ( -AIMAG(BETA2) , REAL ( BET A2) ) 

C = - CEXP(ALENGTH'A4) / CEXP(ALENGTH*DUM1) 

C = C * (BETA 1 + BETAN ' GAMMA * (OMEGA + AMACH * BETA1 )) / ( 

1 BETA2 + BETAN * GAMMA * (OMEGA + AMACH * BETA2 )) 

PS 1 2 = PSI ( DUMMY ) *'2 

AK1=(0MEGA+ AMACH ’BETA1+C* IOMEGA* AMACH*BETA2) ) ’GAMMA 
AK1- CMPLX (-AIMAG(AKI), REAL (AK1J) 

DUM1 = CMPLX ( -A I MAG (BET A 1 ) .REAL (BETA! ) ) 

DUM1 = DUM1 ‘ALENGTH 

DUM2 = CMPLX (-AIMAGIBETA2), REAL (BETA2)) 

DUM2 = DUM2* ALENGTH 

AK2 = ( (OMEGA+AMACH’BETAI ) *CEXP (DUM1 ) + (0MEGA+AMACH*BETA2) *C *CEXP ( 
1DUM2) ) ‘GAMMA 

A<2 = CMPLX (-A1MAG(AK2) .REAL (AK2)) 

CALL CALCA4 ( A4, A5 ) 

BET AW I = (OMEGA* *2 - ETA (LHAT, MHAT, NHAI .ALENGTH, RLMDA) - (AMACH*G2 ( 
INHAT) /EPS I ZN ( LHAT , MHAT , NHAT , ALENGTH , RLMDA ) +BETAN* AK2/EPS I ZN ( LHAT , 

1 MHAT, NHAT, ALENGTH, RLMDA) * (-1 .) ** (NHAT+2) ) /PSI (DUMMY) )/A4 
BETA! = BET AW I - A5 / A4 
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BETA I S 1 1 ) s 8ETAI 

WRITEI6.100) LENGNO 

WRITE (6, 110) LDEX , NOEX 

UR1TEI6, 116) AMACH 

WRITE (6, 1 1 7) GAMMA 

WRITE (6, 1 18) OMEGA 

WRITE(6,119) ALENGTH 

WRITEI6, 101 ) LHAT, MHAT , NHAT 

WRITE (6, 1 15) XI, X2 

WRITE (6, 1 13) F 

WRITE (6. 1 14) A 

WRITEI6, 102) BETAN 

WRITE (6. 120) G 

WRITE (6. 121 ) BETAC 

WRITEI6, 108) AK 

WRITEI6.105) BET AW I 

AN = 1. / GAMMA - BETAWI / AMACH 

WRITE (6, 1 03) AN 

BETAWI = I BETA 1 + C * BETA2 ) / I GAMMA * (OMEGA * ( 1. ♦ C J ♦ 
1 AMACH * ( BETA1 + C * BETA2 ) ) ) 

WRITE (6, 1 12) BETAWI 

AN « 1. / GAMMA - BETAWI / AMACH 

WR!TE(6,105) AN 

WRITE16, 109) BETAISM) 

AN • 1. / GAMMA - BETA IS * 1 ) / AMACH 
WRITEI6, 103) AN 
JX ^ 30 
JY s JX - 1 
CALL CALMU ( JY ) 

JX s JY + 1 

IF ( MATRIX „NE. 3HYES ) GO TO 1 
DO 6 J = 1, JY 
KS = 1 
KT = 3 

WRITE (6, 106) J, KS, KF 
D06 L=1 ,30 

6 WRI TE (6, 1 04) L, ( AMUII.L.J), I = KS, KF ) 

1 CONTINUE 
F = F + FINC 

IF ( JX .EQ. 30 .OR. F .GT. FMAX ) 2 , 7 
50 CALL EXIT 

100 FORMAT ( * 1 THIS IS ENGINE NUMBER *A10) 

101 FORMAT CO THE PRIMARY MODE ASSUMED IS LHAT = M2,* MHAT = \ 12, 



c 

c 


c 

c 


INHAT * M2; 

102 FORMAT ( "0 BET AN * \262U4) 

105 FORMAT ( , +*,70X,* N * •,2621.14." 1*) 

104 FORMAT! 0 * , 13, * M0G13.6) 

105 FORMAT t 0 0 BET AW I s ’,2621. 14) 

106 FORMAT I °0 N J EQUAL M2,' L EQUAL Ml,* TO M2J 

107 FORMAT ( *0° , G21 . 14, 8 WAS THE TIME FOR THE ITERATION *,G21.14) 

108 FORMAT! 0 *’, 70X,° fC s '.2621.14) 

109 FORMAT ( *0 BET A 1 ( 1) - \2G21.14) 

110 FORMAT! 9 **, 70X,* THE COEFFICIENT MATRIX IS ’,12,' X*, 13) 

112 FORMAT 1*0 BETAUID ■ ',2621.14) 

113 FORMAT ! *0 W/WO ■ ’.2G21.14) 

114 FORMAT I '♦* ,70X, ' A * '.2G21.14) 

115 FORMAT ( ’♦’ , 70X, ’ THE LINER BEGINS AT *,F6.4,* AND ENDS AT *,F6.4) 

116 FORMAT ! * 0 THE MACH NUMBER IS *,621.14) 

117 FORMAT !' + ' ,70X, * THE RATIO OF SPECIFIC HEATS IS ',621.14) 

118 FORMAT 1*0 THE FREQUENCY IS ', 2G21.14) 

119 FORMAT!’*', 70X, ' THE LENGTH OF THE COMBUSTOR IS ',621.14) 

120 FORMAT ( '+* , 70X, * G « *, 2621.14) 

121 FORMAT ( * 0 BETAC = ', 2621.14) 

200 FORMAT !A1 0, 7F10.0/6F1 0.0, 513, A5/4F1 0.0) 

END 

FUNCTION BESPRIM (M, L) 

DIMENSION A110, 5) 

**" THESE ARE THE ROOTS OF THE DERIVATIVE OF THE BESSEL FUNCTION OF ORDER ZEO 
**" SET EQUAL TO ZERO 
All ,1) - 0.00000000 
A 12 ,1) = 3.83170597 
A (3 ,1) « 7.01558667 
A 14 ,1) = 10.17546814 
A !5 ,1) - 13.32369194 
A 16 ,1) * 16.47063005 
A (7 ,1) = 19.61585851 
A 18 ,1) * 22.76008438 
A !9 ,1) * 25.90367209 
A! 10,1) « 29.04682853 

"*' THESE are the ROOTS OF THE DERIVATIVE OF THE BESSEL FUNCTION OF ORDER ON 
"" SET EQUAL TO ZERO 
All ,2) = 1.84118378 
A 12 ,2) = 5.33144277 
A (3 ,2) * 8.53631637 
A (4 .2) = 11.70600490 
A (5 ,2) = 14.86358863 
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A (6 ,2) 1 5 ;. 01552786 

A (7 ,2) * 21.16436986 
A (8 ,2) » 24.31132686 
A (9 ,2) * 27.45705057 
A [ 1 0,2) = 30.60192297 

c ,o^ ThESE are ThE roots of the privative OF THE BESSEL FUNCTION OF ORDER Tw 
C 0 9 9 9 SET EQUAL TO ZERO 

All ,3) = 3.05423693 
A (2 ,3) = 6.70613319 
A 13 ,3) = 9.96946782 
A (4 ,3) = 13.17037086 
A 15 ,3) = 16.34752232 
A (6 ,3) = 19.51291278 
A (7 ,3) = 22.67158177 
A (8 ,3) = 25.82603714 
A (9 ,3) * 28.97767277 
A (10,3) = 52.12752702 

C 9 9 9 9 THESE ARE THE ROOTS OF THE DERIVATIVE OF THE BESSEL FUNCTION OF ORDER 3 
C 9999 SET EQUAL TO ZERO 

A(1 ,4) = 4.20118894 
A (2 ,4) = 8.01523660 
A (3 ,4) = 11.34592431 
A (4 ,4) = 14.58584829 
A (5 ,4) = 17.78874787 
A (6 ,4) = 20.97247694 
A (7 ,4) = 24.14489743 
A (8 ,4) = 27.31005793 ' 

A (9 ,4) = 30.47026881 
A ( 1 0 , 4 ) « 33.62694918 

C 9999 THESE ARE THE ROOTS OF THE DERIVATIVE OF THE BESSEL FUNCTION OF ORDER FOR 
C 9999 SET EQUAL TO ZERO 

AM ,5) = 5.31755313 
A (2 ,5) = 9.28239629 
A (3 ,5) = 12.68190844 
A (4 ,5) = 15.96410704 
A (5 ,5) = 19.19602880 
A (6 ,5) = 22.40103227 
A (7 ,5) = 25.58975968 
A (8 ,5) = 28.76783622 
A (9 ,5) = 31.93853934 
A C 1 0 , 5) = 35.10391668 
1 8ESPRIM = A (L, M) 

RETURN 
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ENTRY BESSt'- 

C 9 9 9 9 THESE ARE THE BESSEL NUMBERS Of ORDER ZERO FOR THE ZEROS OF THE BESSEL 

C* 09 * FUNCTION 

AM , 1 ) = 1 . 00000000 
A !2 ,1) = -0.4027588095 
A(3 ,1) = 0.301128303 
A (4 ,1) = -0.249704877 
A (5 ,1) = 0.218359407 
A (6 ,1) = -0.19645371 
A (7 ,1) = 0.180063375 
A (8 ,1) = -0.167184600 
A (9 ,1) = 0.156724985 
A ( 1 0 , 1 ) -0.14801 1108 

C”** THESE ARE THE BESSEL NUMBERS OF ORDER ONE FOR THE ZEROS OF THE BESSEL 

C‘“‘ FUNCTION 

All ,2) = 0.5818649368 
A (2 ,2) = -0.3461258542 
A (3 ,2) = 0.2732981131 
A (4 ,2) = -0.233304416 
A (5 ,2) = 0.207012651 
ACS ,2) = -0.188017488 
A (7 ,2) = 0.173459050 
A 18 ,2) = -0.161838211 
A (9 ,2) = 0.152282069 
AI10, 2) = -0.144242905 

C , ‘’ 9 THESE ARE THE BESSEL NUMBERS OF ORDER TWO FOR THE ZEROS OF THE BESSEL 

C°“ 9 FUNCTION 

Ad ,3) = 0.4864961885 
AI2 ,3) = -0.3135283099 
A (3 ,3) = 0.2547441235 
A (4 ,3) = -0.220881581 
A 15 ,3) = 0.197937434 
A (6 ,3) = -0.161010000 
A (7 ,3) = 0.167835534 
A (8 ,3) = -0.157195167 
A 19 ,3) = 0.148363778 
AI10, 3) = -0.140878333 

THESE ARE THE BESSEL NUMBERS OF ORDER THREE FOR THE ZEROS OF THE BESSEL 
FUNCTION 

Ad ,4) = 0.4343942763 
Af.2 ,4) = -0.2911584413 
A (3 ,4) = 0.240738175 
A (4 ,4) = -0.210965204 
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A 15 ,4) = 0.190419022 
A (6 ,4) = -0.175048405 
A (7 ,4) = 0.162954965 
A (8 ,4) = -0.155102409 
A (9 ,4) = 0.144866574 
A ( 1 0,4) = -0.137844513 

C 00 ** THESE ARE THE BESSEL NUMBERS OE ORDER FOUR FOR THE ZEROS OF THE BESSEL 
C° 9 ** FUNCTION 

All ,5) = 0.3996514545 
A (2 ,5) = -0.2743809949 
A (5 ,5) = 0.229590468 
A (4 ,5) = -0.202763849 
A (5 ,5) = 0.184029896 
A (6 ,5) = -0.169878516 
A (7 ,5) = 0.158655372 
A (8 ,5) - -0.149451156 
A 19 ,5) = 0.141714307 
AI10, 5) = -0.135086328 
GO TO 1 
END 

SUBROUTINE BETAIJU) 

COMPLEX A 1 50, A 151 , A152, AI53, AI55.A155J.DUM1 ,PSI , AMU, BETA IS, COMEG A, 
1DUMMY.FNORM, A, C 1 2 , SUM 1 , SUM2 , SUM3 , SUM4 , AN 
COMPLEX BET A 1 , BETA2, BETAN, BETAC, BETA I , BET AM I , OMEGA, C, AK1 , 

1A<2, PSI2 

COMMON / BLKA / BETA1 , BETA2, AMACH, GAMMA, OMEGA, ALENGTH, C, AK1 
1, AK2, BETAN, RLHDA, BETA], BETAC, BETAMl, LHAT, MHAT, NHAT, 

IPS 12, XI ,X2/BL<B/AMU(1 0,30,29) , BETA IS (30) , A 
COMMON C 1 2 ( 30 , 3C ) , I SET , LDEX , NDEX , ERROR 1 , ERROR2 
NHA1 = NHAT + 1 

COMEGA=CMPLX(-A|MAG (OMEGA) .REAL (OMEGA) ) 

SUM1=SUM2=SUM3=SUM4=CMPLX (0.0.0.0) 

DO 10 LP = 1 , LDEX 
DUM =BESSEL(MHAT+1,LP) 

DO 10 NP = 1 , NDEX 

SUH4=SUM4 + AMU(LP,NP,J-1)*DUM 'C12 (NHAT+1 ,NP) 

10 CONTINUE 

DO 20 NP = 1 , NDEX 
DUM 1 = AMU ( LHAT, NP,J-1 ) 

SUM! = SUM 1 + DUM1 

S'JM2=SUM2 + DUM1 ' (-1 . ) " (NP +1) 

1FINP.EQ.NHAT +1) GOTO 20 

SUM3=SUM3 ♦ DUM1 *2. 'COMEG A* (NP-1 ) ' *2/ 1 (NP-1 ) * *2 -NHAT* *2) 
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I’ll. - (-1 . ) ' * (NP +NHAT -1)) 

20 CONTINUE 

A 1 50 =BET A ] 'AO /PS I (DUMMY ) /EPS I ZN ( LHAT , HHAT , NHAT , ALENGTH, RLMO A ) 

A 151 =BETAN'GAMMA'COMEGA' (M.) "NHAT) *SUM2’FN0RM(LHAT, MHAT, NHAT, 

1 RLMDA , ALENGTH) /EPS I ZN (LHAT , MHAT , NHAT , ALENGTH, RLMD A > 

A 1 52 s (-1 . )' AMACH' SUM3* FNORM (LHAT , MHAT , NHAT , RLMD A , ALENGTH) /EPS 1 ZN ( 

1 LHAT , MHAT , NHAT , ALENGTH, RLMDA) 

A 153=BET AC * FNORM (LHAT , MHAT , NHAT , RLMDA , ALENGTH) 'BESSEL (MHAT+1 .LHAT ) 
1 * SUM4/EPS I RZN (LHAT , MHAT , NHAT , ALENGTH, RLMDA) 

A|55=(AO/PSI (DUMMY) + GAMMA'COMEGA'FNORMILHAT, MHAT, NHAT, RLMDA, 

1 ALENGTH) * SUM 1 ) /EPS I ZN (LHAT , MHAT , NHAT , ALENGTH, RLMDA ) 

BETA1S ( J) = (A 150 -A151 - A152 - A|53)/A155 
WRI TE (6, 1 00) J, BETAIS(J) 

AN = 1. / GAMMA - BETA 1 S ( J) / AMACH 
MR I T£ (6,101) AN 
JS = J 

A 1 50=BET A 1 S ( J) -BETAlSIvM ) 

IF ( ABS (REAL (A 150) /REAL (BETA|S(J-1))).LE. ERROR) .AND. ABS (A1MAG (A 15 
10)/A1MAG(BETA1S(J-1))).LE.ERR0R2) GO TO 2 
RETURN 
2 1 SET = 1 

IF ( REAL (AN) .GT. T . 8 . OR . A I MAG ( AN) .GT. 1 .0JCALL EXIT 
RETURN 

100 FORMAT ( *0 BETA 1 (', 12,*) 8 '.2G21.14) 

101 FORMAT ( ' + ' , 70 X , ’ N = ' ,2621 . 14, ' I*) 

END 

SUBROUTINE CALCA4I A4, A5 ) 

COMPLEX BETA), BETA2, 8ETAN, BETAC, BETA1 , BET AW 1 , OMEGA, C, AO , 

1 AK2, PS12 

COMPLEX A4, A5, PSI, DUMMY, G1, G2 

COMMON / BLKA / BETA), BETA2, AMACH, GAMMA, OMEGA, ALENGTH, C, A<1 
), AK2, BETAN, RLMDA, BET A 1 , BETAC, BET AW I , LHAT, MHAT, NHAT, 

1 PSI2.X1 X2 

A4 = AO/EPS1ZNILHAT, MHAT, NHAT, ALENGTH, RLMDAJ/PS1 (DUMMY) 

A5=BETAC* (G1 ( X2 , NHAT ) -G 1 (X) , NHAT) J/EPSKAT (LHAT, MHAT, NHAT, ALENGTH, 
1RLMDAJ/PS1 (DUMMY) 

RETURN 

END 

SUBROUTINE CALMUI JP ) 

COMPLEX AMU , G 1 , G2 , PS 1 , DUMMY , FNORM , A , DUM3 , C 1 2 
COMPLEX DUM1 , DUM2 , BET A 1 S , SUMS , S , COMEG A , SUM 1 , SUM2 , S’JM3 , SUH4 
COMPLEX BETA), BETA2, BETAN, BETAC, BET A 1 , BET AW I , OMEGA, C, AO, 
)A<2. PS 1 2 
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COMMON / BLKm / BETA! , BET A2 , AMACH, GAMMA, OMEGA, ALENGTH, C, AK1 
1, AK2, BETAN, RLMDA, BETA I, BETAC , BET AW I , LHAT , MHAT, NHAT, 

1PS12, XI .X2/BUB/AMUU 0,30,29) .BETA IS (30) ,A 
COMMON C 1 2 ( 30 , 30 ) , I SET , LDEX , NDEX , ERROR 1 , ERR0R2 
I SET = 0 

DO 50 L = 1 , LDEX 
RLMDA1 = BESPR1M1MHAT+1 ,L) 

AC= BESSEL (MHAT+1 , LHAT) /BESSEL (MHAT+1.LJ 
DO 50 NX = 1 , NDEX 
N=NX- 1 

IFIL.EQ. LHAT. AND. N.EQ.NHAT) GOTO 40 
DUM1 = CMPLX (0. ,0. ) 

IFIL.EQ. LHAT) DUM 1 = (BETA! -BET AW 1 ) 'AK1 / (EPS I ZN (LHAT, MHAT, N, ALENGTH, 
1RLMDA) ) 

DUM1 =A* (BETAC * (G 1 (X2,N)-G1 (XI , N) ) /EPS IKAT (LHAT, MHAT, N, ALENGTH, RLMD 
1 A 1) * AC+DUM1 ) 

DUM2=DUM1/( (OMEGA "2-ETA (LHAT, MHAT, NX, ALENGTH, RLMDAUJ’PS I (DUMMY) ) 
GOTO 50 
40 DUM2=1 . ~ A 

50 AMU (L , NX , 1 ) =DUM2/FN0RM I LHAT , MHAT , NHAT , RLMDA , ALENGTH) 

NX = 1 

COMEGA=CMPLX (-A1MAG (OMEGA) .REAL (OMEGA) ) 

DUM1=CMPLX(0. 0,0.0) 

DUM2=CMPLX (0 . 0,0.0) 

N = NHAT 
NX = NHAT + 1 
DO 2 NY = 1 , NDEX 
NP=NY-1 

IF(N.EQ.NP) 8,7 

7 C12(NX,NY)=GAMMA/2./(N ,, 2-NP ,, 2) * (ALENGTH’COMEGA/3. 1415926' ( (N+NP) 

1 * (S IN ( (N-NP) '3. 1 415926 ’X2/ALENGTH) -$1N( (N-NP) '3. 1 4 1 5926 * X 1 /ALENGTH 
1) ) + (N-NP) ’ (SIN ( (U+NP) *3. 1 4 1 5926 * X2/ALENGTH) -SiNl I N+NP) *3. 1415926‘X 
1 1* ALENGTH) ) )-AMACH*NP* ( (NP+N) * (COS ( (NP-N) ’3. 1 4 1 592S * X2/ALENGTH ) -CO 
IS ( (NP-N) '3. 14 15926'X1 /ALENGTH) )+ (NP-N) ' (COS ( (NP+N) '3. 1 4 1 5926 * X2/AL 
1ENGTH) -COS ( (NP+N) ’3. 1 4 1 5926 * X 1 /ALENGTH) ) ) ) 

GO TO 2 

8 IF (NP.EQ 0)9,10 

S C12(NX,N')=(X2-X1) ‘COMEGA ‘GAMMA 
GO TO 2 

10 C 12 (NX, NY) =GAMMA* (((S1N(2.’N'3. 1 4 1 5926* X2/ALENGTH) -S!N(2. 'N'3. 1415 

1926*X1/ALENGTH))/4./N/3.1415926*ALENGTH+(X2-X1)/2.)’C0MEGA-AMACH/2 
1 . # (SIN (N'3. 1415926*X2/ALENGTH) "2-SIN(N*3. 1415926'XI /ALENGTH) "2) ) 

2 CONTINUE 
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CALL BETA 1 J ( 2 ) 

IF (JP. ECU) RETURN 

00 1 NX 5 1 , NDEX 
N = NX - 1 

DO I NY = 1 , NDEX 
NP=NY-1 

IF(N.EQ.NP) 3,4 

4 C12 (NX , NY) =GAMHA/2./ IN’ *2-NP* ’2) ‘ IALENGTh'COMEGA/3. 1415926 ( (N+NP) 

1 0 (SIN ( (N-NP) *3. 1 4 1 5926 * X2/ALENG TH) -S!N( (N-NP) *3. 1 4 1 5926 * X 1 / ALENGTH 
1 ) )+ (N-NP) * (SIN ( (N+NP) *3. 1415926’X2/ALENGTH) -SINI (N+NP) *3. 1415926*X 
1 1 °ALENGTH) ) ) -AMACH’NP* ( (NP+N) ' (COS ( INP-N) *3. 1 4 1 5926 * X2/ALENGTH) -CO 
IS ( (NP-N) *3. 1 4 1 5926 *X1/ALENGTH) ) + (NP-N) ' (COS ( (NP+N) '3. 1 4 1 5926 * X2/AL 

1 ENG TH ) -C0S( (NP+N) *3. 1 4 1 5926 * X 1 / ALENGTH ) ) ) ) 

GO TO 1 

3 IF (NP.EQ. 0)5,6 

5 C12 (NX,NY) = (X2 - X1 ) 'COMEGA 'GAMMA 
GO TO 1 

6 C 12 (NX, NY) =GAMMA* (( IS IN (2. 'N*3.1415926'X2/ALENGTH) -SIN (2.’N*3. 1415 
1 926 ’X1/ALENGTH) ) /4./N/3. 1 4 1 5926 * ALENG TH+ ( X2-X 1 ) /2. ) 'COMEG A-AMACH/2 
1/ (SIN(N‘3.1415926*X2/ALENGTH)* , 2-SIN(N , 3.1415926'X1/ALENGTH) ,, 2)) 

1 CONTINUE 
DO 60 J=2, JP 
DO 80 LS = 1 , LDEX 
DO 80 NX = 1 , NDEX 
NS=NX- 1 

RLHDA2=BESPRIH( HHAT + 1 , LS) 

DUH3 = BETAC*EPSIZN(LHAT,HHAT, NS , ALENG TH , RLMD A J ‘BESSEL (HHAT+1 , 
1LHAT) /EPS! KATILS, MHAT, NS, ALENGTH.RLMDA2) /BESSEL (MHAT+1 .LS) 

DUM? = CMPLXf 0.0, 0.0 ) 

SUMS=CMPLX(0.,0.) 

IF { LS .EQ. LHAT .AND. NS .EQ. NHAT ) 15, 16 

15 AMU (LHAT, NX, J) = CMPLXI 0.0 , 0.0 ) 

GO TO 80 

16 S'JM1=SUH2=SUM3=SUM4=CMPLX(0. 0,0.0) 

DO 90 LP = 1 , LDEX 
RLMDA3=BESSEL (MHAT +1.LP) 

RLHDA4=BESSEL (MHAT +1 , LHAT ) 

DO 90 NY = 1 , NDEX 

SUN4-SUM4 + AMU(LP,NY, J-1 ) *C 12 (NX, NY) 

50 CONTINUE 

SUM4=SUM4 , DUM3‘RLMDA3/RLMDA4 
DO 21 NP = 1 , NDEX 
DUM1=AM'J(LS,NP,J-1) 
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sum =SUM1 

SUH2=SUH2 ♦ DUH1'M.)"(NP +1) 

IF (NP.EQ.NX) GOTO 21 

SUM3=SUH5 + DUM1 ’2. * C OMEGA * (NP-1 ) "2 * ( 1 .-(-I . ) * ' (NP+NX) ) /( 
1 (NP-1 ) ”2 -NS“2 ) 


21 CONTINUE 

SUM 1 =BET A I S ( J) ‘GAMMA’COMEGA'SUMI 
SUM2=BET AN* GAMMA * COMEG A * (-1 . ) * *NS'SUM2 

SUM3=SUM3 + AMACH’ (3.1415926’NS)”2./2./ALENGTH’AMU(LS,NX,J-1) 
SUM3= 1-1 . ) ' AMACH'SUM3 


DUM1 =CMPLX ( 0 . 0,0.0) 

IF (LS.EQ.LHAT) DUM1= (BETAIStJ) - BETAWI ) * A< 1 
DUM1=DUM1 +DUM3’ (G1 (X2.NS) -GKX1.NS)) 

DUM1 =DUM1 /PS | (DUMMY) /FNORM(LHAT,MHAT,NHAT, RLMDA, ALENGTH) 

AMU (LS, NX, J) = ( DUM 1 +SUM1 +SUM2 +SUM3 ♦ SUM4) / (OMEGA’ *2 - ETA (LHAT 

1 , MHAT.NX, ALENGTH, BESPRiniMHAT+1 ,LSi ) i /'EPSiZNiLHAT.MHAT, NS, ALENGTH, 


2RLMDA) 

80 CONTINUE 

CALL BETA I J ( J + 1 ) 

IF ( I SET .NE. 0 ) GO TO 17 
60 CONTINUE 
17 CONTINUE 
JP s J 
RETURN 
END 

FUNCTION EPS I KATIL.M.N, ALENGTH, RLMDA) 

DUM 1 - ALENGTH 

IF t N .NE. 0) DUM 1 « DUM 1/2. 

IF ( M .EQ. 0) GOTO 1 

EPS KAT = (0.5-M *M / (2. * RLMDA * RLMDA J ) *DUM1 

RETURN 


1 EPS KAT = DUM1/2. 

RETURN 

ENTRY EPS I ZN 

IF ( N .EQ. 0 ) GO TO 2 

EPS1KAT* ALENGTH / 2. 

RETURN 

2 EPS I< AT= ALENGTH 
RETURN 

ENTRY ETA 

EPSIKAT = (N-1 ) ' (N-1 ) ’ 3.1415926 ‘3.1415926 / (ALENGTH ' ALENGTH) 
1+RLMDA • RLMDA 
RETURN 
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end 

FUNCTION EPSIRZN(L,M t N,ALENGTH, RLMDA) 
COMPLEX FN0RMUM2 
DUM1 = 1 . 

IF (M.EQ. 0) DUM 1 = .5 
DUM2=FN0RH(L,M,N,RLMDA,ALENGTH) "2 
EPSIRZN=DUM1 “REAL (DUM2) /3. 1415926 


complex function FNORM (L. m, n, rlmda , alengthj 

DUM 1 = 3.1415926 

IF I H .EQ. 0 ) DUM 1 = 6.2831852 

DUM 1 =DUM 1 “EPSKAT (L,M,N,ALENGTH, RLMDA) ‘BESSEL (M+1 ,L) * '2 
FNORM*CMPLX(DUM1,0.0) 

FNORM* CSQRT (FNORM) 

RETURN 

END 

COMPLEX FUNCTION G1 (X,N) 

COMPLEX BETA! . BETA2, BETAN, BETAC, BETA I, BETAWl , OMEGA, C, AK 1 , 

1 AK2, PS 1 2 

COMPLEX DUM 1 , GDUM, BETA 

COMMON / BLKA / BETA!, BETA2, AMACH, GAMMA, OMEGA, ALE NOTH, C, AK 1 
1, A K2, BETAN, RLMDA, BETA!, BETAC, BETAWl, LHAT, MMAT, NHAT, 

1 PS 12, XI , X2 

GDUM ( A , DUM 1 , BETA . OMEGA , AMACH, X) » ( OMEG A+ AMACH * BET A ) / (A “2-BETA* *2 
1 ) 0 (CEXP(DUMTX) * (DUM 1 *C0S(A*X)+A'S1N(A'X)) -DUM 1 ) 

IF C X .EQ. 0.0 ) GO TO 1 
A = N°3. 1 4 15926/ALENGTH 

BETA = BETA! 

DUM 1 = CMPLX (-AIMAG (BETA! ) .REAL (BETA! ) ) 

G1 = GDUM(A,DUM1, BETA, OMEGA, AMACH, X) 

IF ( N .EQ. 0 .AND. CABS ( BETA2) .LE. 1.E-10 ) GO TO 2 

BETA = BETA2 

DUM! = CMPLX(-AIMAG(BETA2) ,REAL(BETA2)) 

G1 = GAMMA 9 ( G 1 + C * GDUM (A , DUM1 .BETA, OMEGA, AMACH, X) ) 

Cl = CMPLX ( -AIMAG(G1),REAL(GD) 

RETURN 

1 G1 = CMPLX10. 0,0.0) 

RETURN 

2 BETA = C * OMEGA 9 X 

G! = GAMMA 9 ( G1 ♦ BETA ) 

61 = CMPLX ( -AIMAG(GI), REAL (CD) 

RETURN 
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END 

COMPLEX FUNCTION G2(N5 

COMPLEX BETA!, BETA2, BETAN, BETAC, BETA I, BETAWl, OMEGA, C, AK1, 

1 A<2, PS 1 2 

COMPLEX DUM1 , BETA.GBUM, DUM 

COMMON / BLKA / BETA1, BETA2, AMACH, GAMMA, OMEGA, ALENGTH, C, AK1 
1, AK2, BETAN, RLMDA, BETA!, BETAC, BETAWl, LHAT, MHAT, NHAT, 

1 p s 1 2 XI X2 

GBUM(BETA,DUM1,N, AMACH, OMEGA, ALENGTH) = - (AMACH’BETA ’BETA+2 . 'OMEGA 
1 'BETA) / ( ( N ' 3.1415926 / ALENGTH )"2 - BETA "2 ) ' (CEXP (DUMT ' 

1 ALENGTH) ' DUM 1 ' M.)"(N+2) - DUM 1 ) 

DUMIDUM1 ? ALENGTH, NHAT ) = CEXP (D'JMI 'ALENGTH) ' (-1 .)" (NHAT+2) -1 . 

BETA = BETA1 

DUM 1 = CMPLX ( -A ] MAG (BETA 1 ) , REAL (BETA 1 ) ) 

G2 = GBUM (BETA ,DUM1 , N, AMACH, OMEGA, ALENGTH) 

!F( N .EQ. 0 .AND. CABS ( BET A2 ) .LE. 1.E-10 ) RETURN 
BETA = BETA2 

DUM 1 = CMPLX(-AIMAG(BETA2),REAL(BETA2D 

G2 = G2 + C ' GBUM(BETA,DUM1,N, AMACH, OMEGA, ALENGTH) 

RETURN 
ENTRY PSI 

DUM 1 = CMPLX ( -AIMAGIBETA1) , REAL (BET A 1 ) ) 

G2 = DUM 1 'DUMfDUMI .ALENGTH, NHAT! / (NHAT "2' 3. 1415926 "2/ ALENGTH* '2 
1 -BET A 1 * *2) 

IE( NHAT .EQ. 0 .AND. CABS ( BETA2) .LE. 1.E-10 ) GO TO 3 
DUM 1 = CMPLX! -A I MAG (BETA2) , REAL (BETA2) ) 

G2 = ( 62 +C'DUM1'DUM(DUM1, ALENGTH, NHAT) / INHAT"2'3. 1415926**2/ 

1 ALENGTH* '2-BETA2' *2) ) /EPS 1 ZN (LHAT , MHAT , NHAT , ALENGTH , RLMDA) 

RETURN 

3 G2 = ( G2 + C * ALENGTH ) / EPS1ZN(LHAT, MHAT, NHAT, ALENGTH, RLMDA) 
RETURN 
END 


The FOLLOWING IS an EXAMPLE or THE PUNCHED CARD INPUT. 

REFERENCE 0.90 0.00 2.70 0.33 1.20 0.9166 0.00 

5.0 0.0 0.0 0.90 0.01 0.01 1103 30NO 

.0500000 0.0000000 0.9000000 0.0000000 
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THIS IS ENGINE NUMBER REFERENCE 

THE COEFFICIENT MATRIX IS 3X30 

THE MACH NUMBER IS .33000000000000 

THE RATIO OF SPECIFIC HEATS IS 1 . 20000000000C0 

THE FREQUENCY IS 1 .6570654020000 0. 

THE LENGTH OF THE COMBUSTOR IS 2.7000000000000 
THE PRIMARY MODE ASSUMED IS LHAT = 1 MHAT = 1 NHAT 

THE LINER BEGINS AT 0.0000 AND ENDS AT .9000 
W/W0 = .90000000000000 0. 

a = i ; ooooooooooooo o. 

BETAN = 2.74780000000000E-02 0. 

G = .91660000000000 0. 

BETAC = .16666666666667 0. 

K = 5.0000000000000 0. 

BET AW I = .28047988577675 .23813090794998 

N = -1 .66057144750020E-02 -.72160881196963 I 
BET AW I D = .28047988577675 .23813090794997 

N = -1 .660571 44750020E-02 -.72160881196961 1 

BETAM 1) = -3.62072818079238E-02 .16649307640463 

N = .94305236911492 -.50452447395341 I 

BETA 1 ( 2) = 7.574321 55587427E-03 .29698397044242 

N = .81038084437614 -.89995142558308 I 

BETAK 3) = 3.4701 0788505651 E-02 .28016921968467 

N = .72817854893768 -.84899763540808 I 

BETAK 41= 3. 02237853398 196E-02 .27965997230247 

N = .74174610503085 -.8474544615226 2 I 

BETAK 5) = 3.20665167674021E-02 .28192424117869 
N = .73616207040181 -.85431588235965 1 

BETA] ( 6) = 3. 277992041 96444E-02 .28157699072313 

M = .73400024115259 -.83326360825190 I 

BETAK 7) = 3.291 1 69460 35384E- 02 .28169303136132 

N = .73360092544382 -.85361524654944 I 



